GALACTIC DYNAMO

This project is a part of the term project given in course P464 Plamsa Physics and Magnetohydrodynamics taught in Spring 2024 at NISER Bhubaneswar.

Submitted by: Chandan Kumar Sahu, Integrated MSc. SPS batch 19

Supervised by: Dr. Luke R. Chamandy, SPS, NISER

---

---

CHAPTER 1

Contents

1. Problem Statement
2. Theory
   1. The mean-field induction equation
   2. Solving in cylindrical coordinates
   3. Calculation of total magnetic field magnitude and pitch angle
   4. Magnetic decay constant
3. Numerical solution to Galactic Magnetic Fields
   1. Solver: Crank-Nicholson
   2. Defining the grid
   3. Boundary conditions

THE CODE

1. Example 1 - Simplest case ($B_r>0$ and $B_\phi>0$ at all times)
2. Example 2 - ($B_\phi>0$ at all times)
3. Example 3 - Arbitrary seed field

---

---

1.1 - Problem Statement

Solve the diffusion equation in the z-direction.

• Explore the evolution of the magnetic field magnitude and of the exponential decay rate.
• Explore the evolution of the spatial solution for $B_r$ and $B_\phi$, and of the pitch angle of the mean magnetic field p.
• Explore how different boundary conditions affect the results.
• Explore how different seed fields affect the results.

---

---

1.2 - Theory

The dynamo theory describes the process through which a rotating, convecting, and electrically conducting fluid can maintain a magnetic field over astronomical time scales. The evolution of the magnetic field is described by the induction equation, which relates changes in the magnetic field to the velocity field of the conducting fluid and its magnetic diffusivity.

In this chapter, we will describe how the galactic magnetic field varies with time in the $z$-direction, and discuss the points mentioned in the problem statement.

1.2.1 - The mean-field induction equation

We have the mean-field induction equation as $$ \dfrac{\partial \bar{\mathbf{B}}}{\partial t} = \nabla \times \left[ \bar{\mathbf{V}} \times \bar{\mathbf{B}} + \mathcal{E} - \eta \left( \nabla \times \bar{\mathbf{B}} \right) \right] $$ where $\mathcal{E} = \left( \alpha \bar{\mathbf{B}} \right) - \eta_t \left( \nabla \times \bar{\mathbf{B}} \right)$

We will solve the equations in the cylindrical coordinates (r, $\phi$, z) with the origin at the galactic centre and the z-axis parallel to the galactic angular velocity. However, to simplify things, lets make some approximations.

1. Omit the terms involving $\bar{\mathbf{V}} \times \bar{\mathbf{B}}$ and $\alpha$. We will land up on just the diffusion equation.

$$ \dfrac{\partial \bar{\mathbf{B}}}{\partial t} = - \nabla \times \left[ \eta_T \left( \nabla \times \bar{\mathbf{B}} \right) \right] $$

where $\eta_T = \eta + \eta_t$

2. Take $\eta_T$ independent of $\bar{\mathbf{B}}$. Our equation becomes

$$ \dfrac{\partial \bar{\mathbf{B}}}{\partial t} = - \eta_T \left[ \nabla \times \left( \nabla \times \bar{\mathbf{B}} \right) \right] $$

But $\nabla \times \left( \nabla \times \bar{\mathbf{B}} \right) = \nabla \left( \nabla \cdot \bar{\mathbf{B}} \right) - \nabla^2 \bar{\mathbf{B}} $ and $\nabla \cdot \bar{\mathbf{B}} = 0$ (Gauss's Law), so we finally have $$ \boxed{ \dfrac{\partial \bar{\mathbf{B}}}{\partial t} = \eta_T \nabla^2 \bar{\mathbf{B}} }$$

This is the Fickian diffusion equation. We will solve this equation numerically.

---

1.2.2 - Solving in cylindrical coordinates

In cylindrical coordinates, $$ \mathbf{\bar{B}} = \bar{B}_r \mathbf{\hat{r}} + \bar{B}_{\phi} \mathbf{\hat{\phi}} + \bar{B}_z \mathbf{\hat{z}} $$

The $\nabla^2B$ becomes $$ \begin{aligned} \nabla^2 \mathbf{\bar{B}} = & \left[ \frac{\partial}{\partial r} \left( \frac{1}{r} \frac{\partial}{\partial r} \left( r \bar{B}_r \right) \right) + \frac{1}{r^2} \frac{\partial^2 \bar{B}_r}{\partial \phi^2} + \frac{\partial^2 \bar{B}_r}{\partial z^2} - \frac{2}{r^2} \frac{\partial \bar{B}_\phi}{\partial \phi} \right] \mathbf{\hat{r}} \\ & + \left[ \frac{\partial}{\partial r} \left( \frac{1}{r} \frac{\partial}{\partial r} \left( r \bar{B}_\phi \right) \right) + \frac{1}{r^2} \frac{\partial^2 \bar{B}_\phi}{\partial \phi^2} + \frac{\partial^2 \bar{B}_\phi}{\partial z^2}+\frac{2}{r^2} \frac{\partial \bar{B}_r}{\partial \phi} \right] \mathbf{\hat{\phi}} \\ & + \left[\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial \bar{B}_z}{\partial r}\right)+\frac{1}{r^2} \frac{\partial^2 \bar{B}_z}{\partial \phi^2}+\frac{\partial^2 \bar{B}_z}{\partial z^2}\right] \mathbf{\hat{z}} \end{aligned} $$

So we get the final component-wise equations as $$ \begin{aligned} \frac{\partial \bar{B}_r}{\partial t} &= \eta_T \left[ \frac{\partial}{\partial r} \left( \frac{1}{r} \frac{\partial}{\partial r} \left( r \bar{B}_r \right) \right) + \frac{1}{r^2} \frac{\partial^2 \bar{B}_r}{\partial \phi^2} + \frac{\partial^2 \bar{B}_r}{\partial z^2} - \frac{2}{r^2} \frac{\partial \bar{B}_\phi}{\partial \phi} \right] \\ \frac{\partial \bar{B}_\phi}{\partial t} &= \eta_T \left[ \frac{\partial}{\partial r} \left( \frac{1}{r} \frac{\partial}{\partial r} \left( r \bar{B}_\phi \right) \right) + \frac{1}{r^2} \frac{\partial^2 \bar{B}_\phi}{\partial \phi^2} + \frac{\partial^2 \bar{B}_\phi}{\partial z^2}+\frac{2}{r^2} \frac{\partial \bar{B}_r}{\partial \phi} \right] \\ \frac{\partial \bar{B}_z}{\partial t} &= \eta_T \left[\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial \bar{B}_z}{\partial r}\right)+\frac{1}{r^2} \frac{\partial^2 \bar{B}_z}{\partial \phi^2}+\frac{\partial^2 \bar{B}_z}{\partial z^2}\right] \end{aligned} $$

Since the problem statement in the project aims to solve only in the z-direction, we remove all radial or azimuthal variations of the magnetic field $\left(\dfrac{\partial }{\partial r} = \dfrac{\partial }{\partial \phi} = 0 \right)$.

We are now left with these simple equations to solve, i.e., the Fickian diffusion equations. $$ \boxed{ \frac{\partial \bar{B}_r}{\partial t} = \eta_T \frac{\partial^2 \bar{B}_r}{\partial z^2} } \qquad \qquad \qquad \boxed{ \frac{\partial \bar{B}_\phi}{\partial t} = \eta_T \frac{\partial^2 \bar{B}_\phi}{\partial z^2} } \qquad \qquad \qquad \boxed{ \frac{\partial \bar{B}_z}{\partial t} = \eta_T \frac{\partial^2 \bar{B}_z}{\partial z^2} } $$

---

1.2.3 - Calculation of total magnetic field magnitude ($B_{total}$) and pitch angle ($p_B$)

We can calculate the magnitude of the total magnetic field as $$ B_{\text{total}} = \sqrt{\bar{B}_r^2 + \bar{B}_\phi^2} $$ And $$ p_B = \tan^{-1} \left( \dfrac{\bar{B}_r}{\bar{B}_\phi} \right) \qquad \qquad \text{where} \qquad -\dfrac{\pi}{2} < p_B < \dfrac{\pi}{2}$$

---

1.2.4 - Magnetic decay constant

The diffusion of the total magnetic field can be expressed in the form

$$ B_{\text{total}}(z, t) = \tilde{B}(r)\exp(\gamma t) $$

where $\tilde{B}(r)$ contains all the variation in $r$ and the exponential factor contains time variation, where $\gamma$ is the magnetic decay constant.

Our final goal is to calculate this $\gamma$.

---

---

1.3 - Numerical Solution to Galactic magnetic Fields

Now we begin to solve the diffusion equations numerically.

We will first solve the heat diffusion equation, and use them to calculate the total $B_{\text{total}}$ and the pitch angle $p_B$. Finally we will calculate the decay constant $\gamma$.

1.3.1 - Solver: Crank-Nicholson

We use the Crank-Nicholson algorithm to solve the diffusion equation. It combines implicit and explicit schemes, resulting in a stable and accurate solution.

Let us now get into the details of discretization. Given below is the discretization scheme for the forward Euler explicit method.

The diffusion equation is given as $$ \frac{\partial B}{\partial t} = \eta_T \frac{\partial^2 B}{\partial z^2} $$ We use numerical derivative scheme for discretization. The superscript $n$ denotes index in time while the subscript $i$ denotes index in space. $dt$ and $dz$ denote cell thickness in time and space respectively. $$ \frac{B^{j+1}_{i} - B^{j}_{i}}{dt} = \eta_T \: \frac{B^{j}_{i+1} - 2B^{j}_{i} + B^{j}_{i-1}}{dz^2} $$ This solves to $$ B^{j+1}_{i} = B^{j}_{i} + \dfrac{\eta_T dt}{dz^2} \: \left(B^{j}_{i+1} - 2B^{j}_{i} + B^{j}_{i-1} \right) $$ which gives the solution B at the next time-step using the information from the previous time-step. A stencil of forward Euler method is shown below:

Forward Euler Method

Crank Nicholson Method

The stability of explicit Euler methods require $\dfrac{\eta_T dt}{dz^2} < 0.5$ which has the disadvantages of instability or high computation power. To overcome this, implicit schemes were developed, which use the information from both previous and present time-step to obtain values at the repsent time-step.

In this project, we will use the Crank-Nicholson method which uses average of previous and present time-step to solve the problem. Here is the discretization scheme for it.

$$ \frac{B^{j+1}_{i} - B^{j}_{i}}{dt} = \dfrac{\eta_T}{2} \: \left( \dfrac{B^{j+1}_{i+1} - 2B^{j+1}_{i} + B^{j+1}_{i-1}}{dz^2} \right) + \dfrac{\eta_T}{2} \: \left( \dfrac{B^{j}_{i+1} - 2B^{j}_{i} + B^{j}_{i-1}}{dz^2} \right) $$

Now we separate the present time-step $(j+1)$ and the past time-step $(j)$ $$ B^{j+1}_{i} - \dfrac{\eta_T \: dt}{2 \: dz^2} \: \left( B^{j+1}_{i+1} - 2B^{j+1}_{i} + B^{j+1}_{i-1} \right) = B^{j}_{i} + \dfrac{\eta_T \: dt}{2 \: dz^2} \: \left( B^{j}_{i+1} - 2B^{j}_{i} + B^{j}_{i-1} \right) $$

Taking $\dfrac{\eta_T \: dt}{2 \: dz^2} = \sigma$, we rewrite the equations as $$ \left(1+2\sigma \right) B^{j+1}_{i} - \sigma B^{j+1}_{i+1} - \sigma B^{j+1}_{i-1} = \left(1-2\sigma \right) B^{j}_{i} + \sigma B^{j}_{i+1} + \sigma B^{j}_{i-1} $$

We can now form a matrix for the two timesteps and solve it. $$ MB^{j+1} = NB^{j} $$ where $$ M=\left[\begin{array}{ccccc}1+2 \sigma & -\sigma & 0 & \cdots & 0 \\ -\sigma & 1+2 \sigma & -\sigma & \ddots & \vdots \\ 0 & -\sigma & 1+2 \sigma & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & -\sigma \\ 0 & \cdots & 0 & -\sigma & 1+2 \sigma\end{array}\right] \qquad \qquad N=\left[\begin{array}{ccccc}1-2 \sigma & \sigma & 0 & \cdots & 0 \\ \sigma & 1-2 \sigma & \sigma & \ddots & \vdots \\ 0 & \sigma & 1-2 \sigma & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \sigma \\ 0 & \cdots & 0 & \sigma & 1-2 \sigma\end{array}\right] $$

The final answer i.e., $B^{j+1}$ is obtained as $$ B^{j+1} = M^{-1}NB^{j} $$

By averaging the values of variables at current and next time steps, it reduces numerical errors and suppresses oscillations.

The method is highly stable, allowing for larger time steps compared to explicit methods, and is second-order accurate in time and space. This makes it well-suited for simulating heat diffusion processes with high accuracy and efficiency.

---

1.3.2 - Defining the grid

$\star$ Spatial Grid ($z$)

The size of a galactic disk typically ranges from around 10 kpc to 50 kpc. Our Milky way has a radius of ~40 kpc. Within this disk, there are thin and thick disks. The thin disk, where most of the younger stars reside, typically extends from the galactic center to about 100-1000 pc. Meanwhile, the thick disk, containing older stars and having a higher vertical velocity dispersion, can extend further out to ~3 kpc.

Considering these typical values, we choose the thickness of our model galaxy to be ~200 pc.

So, $z$ extends from -100 pc to +100 pc. For solving numerically, we normalize the spatial grid to -1 to 1, with a grid cell size of $dz = 0.01$.

Figure: Schematic of the Milky Way galaxy

$\star$ Temporal Grid ($t$)

The typical variation of magnetic fields usually vary in millions of years. However, the duration of simulation depends on the seed field, so we change this values and its binning according the the seed field. In general, we take 100 points within the temporal grid range.

$\star$ Magnetic diffusion coefficient ($\eta_T$)

The magnetic diffusion coefficient measures the rate at which magnetic fields diffuse through a medium. We calculate the value of $\eta_T$ as $$ \eta_T \approx \dfrac{1}{3} \tau v_{\text{rms}}^2 $$

But we first need to calculate its value in the normlizated units. For a typical galaxy, the faraday time $\tau \approx 10$ Myr and velocity $v_{\text{rms}} \approx 10^4$ km/s

$$ \begin{aligned} \eta_T &= \dfrac{1}{3} \times 10 \text{ Myrs} \times (10 \text{ km/s})^2 \\ &= \dfrac{100}{3} \times \text{ Myrs} \times \dfrac{\text{ km}^2}{\text{s}^2} \\ &= \dfrac{100}{3} \times \text{ Myrs} \times \left( \dfrac{\text{ km} \times 100 \text{ pc}}{3.086\times 10^{15}\text{ km}} \right)^2 \left(\dfrac{3.154\times 10^{13} \text{ s}}{\text{s} \times \text{Myr}} \right)^2 \\ &= \dfrac{100}{3} \times \left( \dfrac{3.154\times 10^{13}}{3.086\times 10^{15}} \right)^2 \times \text{ (100 pc)}^2/\text{ Myr} \\ &= 3.48 \times 10^{-3}\text{ (100 pc)}^2/\text{ Myr} \end{aligned} $$

Hence, we use the following values $$ -1<z<1 \qquad \qquad \text{ with } dz=0.01 \text{ and } \qquad \qquad \eta_T = 3.48 \times 10^{-3}\text{ (100 pc)}^2/\text{ Myr} $$ Here we have $h = 100$ pc

---

1.3.3 - Boundary conditions

Two common types of boundary conditions are vacuum boundary conditions and isolated boundary conditions.

Vacuum Boundary Conditions (BC):

In vacuum boundary conditions, the system is assumed to be in contact with a perfect vacuum, the boundaries of the system under consideration are devoid of any matter or energy. It ensures that the quantity becomes 0 at the boundaries at all time.

Isolated Boundary Conditions:

Isolated boundary conditions imply that the system is completely isolated from its surroundings, meaning there are no exchanges of matter or energy across its boundaries. It ensures that the derivative of the quantity becomes 0 at the boundaries at all time.

We are calculating the diffusion of $B_r$ and $B_\phi$. As the distance in z-increases beyond, $B_z$ dominates over $B_r$.

Within the thin disk approximation, it is safe to assume that $B_r \rightarrow 0$ as $z \rightarrow 100$ pc. The same thing applies to $B_\phi$ too. Hence, we apply vacuum boundary conditions for all the cases here.

---

---

THE CODE

For the full code with al functions and plotting scripts, please visit my github repository.

Now let us get into some examples to understand how the magnetic field varies with different types of seed fields.

from my_code import *
from plotting import *

Example 1

Let us begin with simple seed fields, something which has a central peak and vanishes at the boundary.

Defining the seed field

def init_cond_Br(x):
    return 1-x**2

def init_cond_Bphi(x):
    return np.cos(np.pi/2*x)**2

def source_term(x, t):
    return 0

z = np.linspace(-1, 1, 101)
title_1 = r'Parabolic seed field '+'\n'+r'$ y = 1-x^2$'
title_2 = r'Cos square seed field '+'\n'+r'$ y = \cos^2(\pi x/2)$'

plot_init_cond(z, init_cond_Br, init_cond_Bphi, title_1, title_2, global_title='INITIAL CONDITIONS FOR MAGNETIC FIELD COMPONENTS')
plt.show()

Solving the diffusion equation

# Constants and parameters
eta_T = 3.48e-3    # magnetic diffusivity
t_max = 200     # total simulation time
z_min = -1.0     # minimum thickness of the disc
z_max = 1.0     # thickness of the disc
dt = t_max/200       # time step
dz = 0.01       # spatial step in z direction

# Solve the diffusion equation in radial direction
solution_r, spatial_grid, time_grid = crank_nicolson_diffusion(z_min, z_max, t_max, dz, dt, eta_T, init_cond_Br, source_term, diff_matrix_vacuum_boundary)
# Solve the diffusion equation in azimuthal direction
solution_phi, spatial_grid, time_grid = crank_nicolson_diffusion(z_min, z_max, t_max, dz, dt, eta_T, init_cond_Bphi, source_term, diff_matrix_vacuum_boundary)

# Plot the diffusion equation solution
plot_diff(time_grid, spatial_grid, solution_r, solution_phi)
plt.show()

# print(time_grid)

Calculation of $B_{\text{total}}$ and $p_B$

B_total, pitch = get_B_and_pitch(solution_r, solution_phi)

# Plot the total magnetic field and the pitch angle
plot_pitch(time_grid, spatial_grid, B_total, pitch)
plt.show()

Decay rate calculation

# Plot the decay of the magnetic field at the midplane
B_mid = np.log(B_total[int(len(spatial_grid)/2), :])
m, c = np.polyfit(time_grid[-50:], B_mid[-50:], 1)

plot_decay(time_grid, B_mid, m, c)
plt.show()

So the value of magnetic decay rate or decay constant obtained for this case is $\gamma = -8.417 \times 10^{-3}$

Example 2

Now we will move to slightly complex seed fields, where the variations and could be arbitrary in both $r$ and $\phi$. However, we have chosen the seed $B_\phi$ to be always positive in this exmaple.

Negative magnitudes correspond to reversal in the direction of the magnetic field.

Defining the seed field

def init_cond_Br(x):
    return -2*x*np.sin(2*np.pi*x)

def init_cond_Bphi(x):
    return 0.75 - 0.5*x**2 + 0.25*np.cos(3*np.pi*x)

def source_term(x, t):
    return 0

z = np.linspace(-1, 1, 101)
title_1 = r'Algebraic-sinusoidal seed field '+'\n'+r'$ y = -2x\sin(2\pi x)$'
title_2 = r'Sinusoidal seed field '+'\n'+r'$ y = \sin(2\pi x)$'

plot_init_cond(z, init_cond_Br, init_cond_Bphi, title_1, title_2, global_title='INITIAL CONDITIONS FOR MAGNETIC FIELD COMPONENTS')
plt.show()

# Constants and parameters
eta_T = 3.48e-3    # magnetic diffusivity
t_max = 200     # total simulation time
z_min = -1.0     # minimum thickness of the disc
z_max = 1.0     # thickness of the disc
dt = t_max/200       # time step
dz = 0.01       # spatial step in z direction

# Solve the diffusion equation in radial direction
solution_r, spatial_grid, time_grid = crank_nicolson_diffusion(z_min, z_max, t_max, dz, dt, eta_T, init_cond_Br, source_term, diff_matrix_vacuum_boundary)
# Solve the diffusion equation in azimuthal direction
solution_phi, spatial_grid, time_grid = crank_nicolson_diffusion(z_min, z_max, t_max, dz, dt, eta_T, init_cond_Bphi, source_term, diff_matrix_vacuum_boundary)

# Plot the diffusion equation solution
plot_diff(time_grid, spatial_grid, solution_r, solution_phi)

plt.show()

Calculation of $B_{\text{total}}$ and $p_B$

B_total, pitch = get_B_and_pitch(solution_r, solution_phi)

# Plot the total magnetic field and the pitch angle
plot_pitch(time_grid, spatial_grid, B_total, pitch)
plt.show()

Decay rate calculation

# Plot the decay of the magnetic field at the midplane
B_mid = np.log(B_total[int(len(spatial_grid)/2), :])
m, c = np.polyfit(time_grid[-50:], B_mid[-50:], 1)

plot_decay(time_grid, B_mid, m, c)
plt.show()

So the value of magnetic decay rate or decay constant obtained for this case is $\gamma = -8.417 \times 10^{-3}$

Example 3

Now let us move to completely arbitrary seed field spanning both positive and negative values.

Defining the seed field

def init_cond_Br(x):
    # return 0.75 - 0.5*x**2 + 0.25*np.cos(5*np.pi*x)
    return (1-x**2)*np.cos(2*np.pi*x)

def init_cond_Bphi(x):
    return (1-x**2)*np.cos(np.pi*x)

def source_term(x, t):
    return 0

z = np.linspace(-1, 1, 101)
title_1 = r'Arbitrary seed field '+'\n'+r'$ y = 0.75 - 0.5x^2 + 0.25\cos(5\pi x)$'
title_2 = r'Algebraic-sinusoidal seed field '+'\n'+r'$ y = -2x\sin(2\pi x)$'

plot_init_cond(z, init_cond_Br, init_cond_Bphi, title_1, title_2, global_title='INITIAL CONDITIONS FOR MAGNETIC FIELD COMPONENTS')
plt.show()

# Constants and parameters
eta_T = 3.48e-3    # magnetic diffusivity
t_max = 200     # total simulation time
z_min = -1.0     # minimum thickness of the disc
z_max = 1.0     # thickness of the disc
dt = t_max/200       # time step
dz = 0.01       # spatial step in z direction

# Solve the diffusion equation in radial direction
solution_r, spatial_grid, time_grid = crank_nicolson_diffusion(z_min, z_max, t_max, dz, dt, eta_T, init_cond_Br, source_term, diff_matrix_vacuum_boundary)
# Solve the diffusion equation in azimuthal direction
solution_phi, spatial_grid, time_grid = crank_nicolson_diffusion(z_min, z_max, t_max, dz, dt, eta_T, init_cond_Bphi, source_term, diff_matrix_vacuum_boundary)

# Plot the diffusion equation solution
plot_diff(time_grid, spatial_grid, solution_r, solution_phi)

plt.show()

Calculation of magnitude of $B_{total}$ and $p_B$

B_total, pitch = get_B_and_pitch(solution_r, solution_phi)

# Plot the total magnetic field and the pitch angle
plot_pitch(time_grid, spatial_grid, B_total, pitch)
plt.show()

One observes that there are the variations in the pitch angle are abrupt in this case. This is because there is change in the direction of the magnetic fields, causing the pitch angle to move from $+90^\circ$ to $-90^\circ$, which is natural. One must not think of this as code instability.

Decay rate calculation

# Plot the decay of the magnetic field at the midplane
B_mid = np.log(B_total[int(len(spatial_grid)/2), :])
m, c = np.polyfit(time_grid[-50:], B_mid[-50:], 1)

plot_decay(time_grid, B_mid, m, c)
plt.show()

So the value of magnetic decay rate or decay constant obtained for this case is $\gamma = -8.418 \times 10^{-3}$

---

---

CHAPTER 2

Contents

1. Problem Statement
2. Theory
   1. Revisiting the mean-field induction equation
   2. Dynamo number
3. Numerical Solution to Galactic magnetic Fields with $\alpha-\Omega$ approximation
   1. Solver: Crank-Nicholson
   2. Grid and boundary conditions

THE CODE

1. Example 1 - Decaying fields
2. Example 2 - Growing fields

---

---

2.1 - Problem Statement

Solve the mean-field $\alpha-\Omega$ dynamo equations in the kinematic regime. That is, include the $\Omega$ effect term in the equation for $\dfrac{\partial \bar{B}_\phi}{\partial t}$ and the $\alpha$ effect term in the equation for $\dfrac{\partial \bar{B}_r}{\partial t}$. This requires specifying the overall magnitude and spatial dependence of $\Omega$ and $\alpha$.

• Repeat the investigation you had done for task 1, with the new equations, for different values of the dynamo number, $D$

$$ D = − \dfrac{\alpha_0 q \Omega h^3}{\eta_t^2} $$

where $q = − \dfrac{d \ln \Omega}{d \ln r}$ and $\alpha_0 > 0$ is the amplitude of the $\alpha$ effect. Note that $q > 0$ if $\Omega$ decreases with $r$, which is generally the case in galaxies, so $D < 0$.

• The exponential decay becomes exponential growth if $|D| > |D_c|$, where $D_c$ is the critical dynamo number. Find the critical dynamo number numerically (ideally, you would automate this feature).
• Compare the growth rate you obtain for a given value of $D$ (for $|D| > |D_c|$) with the no-z solution prediction for the local growth rate $\gamma$. Do the same comparison for $D_c$ . Do the results agree with your expectations?

---

---

2.2 - Theory

In the previous chapter, we studied how the galactic magnetic field varies with time in the $z$-direction. For that we covered the diffusion equation, its numerical solution using Crank-Nicholson method and some example seed fields, and finally calculated the decay rate.

2.2.1 - Revisiting the mean-field induction equation

In this chapter, we will begin with the kinematic regime of the solutions, beginning with the addition of the $\alpha-\Omega$ term in uor equations. Lets revisit the mean-field induction equation.

$$ \dfrac{\partial \bar{\mathbf{B}}}{\partial t} = \nabla \times \left[ \bar{\mathbf{V}} \times \bar{\mathbf{B}} + \mathcal{E} - \eta \left( \nabla \times \bar{\mathbf{B}} \right) \right] $$

where $\mathcal{E} = \left( \alpha \bar{\mathbf{B}} \right) - \eta_t \left( \nabla \times \bar{\mathbf{B}} \right)$

We will solve the equations in the cylindrical coordinates (r, $\phi$, z) with the origin at the galactic centre and the z-axis parallel to the galactic angular velocity. However, to simplify things, lets make some approximations again.

1. Consider $\bar{V}$ to be of this form, taking $\eta_T = \eta + \eta_t$,

$$ \mathbf{\bar{V}} = \bar{V}_r(r) \mathbf{\hat{r}} + r \Omega(r) \mathbf{\hat{\phi}} + \bar{V}_z(r, z) \mathbf{\hat{z}} $$

2. Take $\eta_T$ independent of $\bar{\mathbf{B}}$. Our equation becomes

But $\nabla \times \left( \nabla \times \bar{\mathbf{B}} \right) = \nabla \left( \nabla \cdot \bar{\mathbf{B}} \right) - \nabla^2 \bar{\mathbf{B}} $ and $\nabla \cdot \bar{\mathbf{B}} = 0$ (Gauss's Law). So we get,

$$ \boxed{ \dfrac{\partial \bar{\mathbf{B}}}{\partial t} = \nabla \times \left( \bar{\mathbf{V}} \times \bar{\mathbf{B}} \right) + \nabla \times \left(\alpha \bar{\mathbf{B}} \right) - \eta_T \left( \nabla \times \nabla \times \bar{\mathbf{B}} \right) } $$

3. Solving again in cylindrical coordinates and taking azimuthal symmetry, we can write the equations for $\bar{B}_r$ and $\bar{B}_{\phi}$ as

$$ \frac{\partial \bar{B}_r}{\partial t} = V_r \frac{\partial \bar{B}_z}{\partial z} - \frac{\partial (V_z \bar{B}_r)}{\partial z} - \frac{\partial (\alpha \bar{B}_\phi)}{\partial z} + \eta_T \left[ \frac{\partial^2 \bar{B}_r}{\partial z^2} + \frac{\partial}{\partial r} \left( \frac{1}{r} \frac{\partial}{\partial r} \left( r \bar{B}_r \right) \right) \right] $$$$ \frac{\partial \bar{B}_\phi}{\partial t} = r \Omega \frac{\partial \bar{B}_z}{\partial z} - \frac{\partial (V_z \bar{B}_\phi)}{\partial z} - \frac{\partial (V_r \bar{B}_\phi)}{\partial r} + \frac{\partial (r \Omega \bar{B}_r)}{\partial r} + \frac{\partial (\alpha \bar{B}_r)}{\partial z} - \frac{\partial (\alpha \bar{B}_z)}{\partial r} + \eta_T \left[ \frac{\partial^2 \bar{B}_\phi}{\partial z^2} + \frac{\partial}{\partial r} \left( \frac{1}{r} \frac{\partial}{\partial r} \left( r \bar{B}_\phi \right) \right) \right] $$

4. Taking the $\alpha-\Omega$ approximation and putting the value $q = -\dfrac{\partial \ln \Omega}{\partial \ln r} = -\dfrac{r}{\Omega} \dfrac{\partial \Omega}{\partial r}$ we obtain

$$ \frac{\partial \bar{B}_r}{\partial t} = - \frac{\partial (\alpha \bar{B}_\phi)}{\partial z} + \eta_T \frac{\partial^2 \bar{B}_r}{\partial z^2} $$$$ \frac{\partial \bar{B}_\phi}{\partial t} = -q \Omega \bar{B}_r + \eta_T \frac{\partial^2 \bar{B}_\phi}{\partial z^2} $$

Here we have neglected all terms which involved $\dfrac{\partial}{\partial r}$ in $B$. We have also omitted the $\alpha^2$ term for now.

5. Finally consider $\alpha = \alpha_0$ to be independent of $z$. So our final equations are

$$ \boxed{ \frac{\partial \bar{B}_r}{\partial t} = - \alpha_0\frac{\partial (\bar{B}_\phi)}{\partial z} + \eta_T \frac{\partial^2 \bar{B}_r}{\partial z^2}} \qquad \qquad \text{and} \qquad \qquad \boxed{ \frac{\partial \bar{B}_\phi}{\partial t} = -q \Omega \bar{B}_r + \eta_T \frac{\partial^2 \bar{B}_\phi}{\partial z^2} } $$

---

2.2.2 - Dynamo number

The dynamo number is defined as

$$ D = − \dfrac{\alpha_0 q \Omega h^3}{\eta_t^2} $$

Here $h$ is the thickness of the thin disk of the galaxy, $\alpha$ is the term responsible for the twisting of the toroidal fields into poloidal fields, $\Omega$ is the rotation rate fo the galaxy, causing the twisting of poloidal fields into toroidal fields, and $q = -\dfrac{\partial \ln \Omega}{\partial \ln r} = -\dfrac{r}{\Omega} \dfrac{\partial \Omega}{\partial r}$.

One can parameterize $\Omega$ as $$ \Omega = \dfrac{\Omega_0}{\sqrt{1 + \left( \dfrac{r}{r_0} \right)^2 }} $$ where $r_0$ is the radius of the galaxy.

Thus $q$ solves as $$ q = - \left( \dfrac{r^2}{r_0^2} \right) \left( 1 + \dfrac{r^2}{r_0^2} \right)^{-1} $$

With the values $\eta_T = 3.48 \times 10^{-3}\text{ (100 pc)}^2/\text{ Myr}$ (as chosen previously),

$r_0 = 30$ kpc (typical radius of a galaxy),

$h = 100$ pc (typical thickness of thin disks),

$\Omega = 40$ km/s/kpc (typical order of magnitude rotation rate of galaxies)

we get $q = -0.5$

From this, the critical dynamo number is found to be, $$ \begin{aligned} D &= − \dfrac{\alpha_0 \times -0.5 \times 40 \text{ km/s/kpc} \times (100 \text{ pc})^3}{(3.48 \times 10^{-3}\text{ (100 pc)}^2/\text{ Myr})^2} \\ &= 1.725 \times 10^7 \: \alpha_0 \end{aligned} $$

where $\alpha_0$ is in km/s.

2.3 - Numerical Solution to Galactic magnetic Fields with $\alpha-\Omega$ approximation

Now we begin to solve the equations numerically.

2.3.1 - Solver: Crank-Nicholson (modified for coupled equations)

We had used the Crank-Nicholson method for solving the 1D diffusion equation in the previous chapter, however, things get complicated when we mpove to coupled diffusion equations, like in this case. Let us understand the discretization scheme for coupled diffusion equations.

Let $\bar{B}_r = P$ and $\bar{B}_\phi = Q$ for understanding purposes. Our equations modify to $$ \frac{\partial P}{\partial t} = - \alpha_0 \frac{\partial Q}{\partial z} + \eta_T \frac{\partial^2 P}{\partial z^2} $$ $$ \frac{\partial Q}{\partial t} = -q \Omega P + \eta_T \frac{\partial^2 Q}{\partial z^2} $$

Lets first discretize equations, $$ \frac{P^{j+1}_{i} - P^{j}_{i}}{dt} = -\dfrac{\alpha_0}{2} \: \left( \frac{Q^{j+1}_{i+1} - Q^{j+1}_{i}}{dz} + \frac{Q^{j}_{i+1} - Q^{j}_{i}}{dz} \right) \dfrac{\eta_T}{2} \: \left( \dfrac{P^{j+1}_{i+1} - 2P^{j+1}_{i} + P^{j+1}_{i-1}}{dz^2} + \dfrac{P^{j}_{i+1} - 2P^{j}_{i} + P^{j}_{i-1}}{dz^2} \right) $$ and $$ \frac{Q^{j+1}_{i} - Q^{j}_{i}}{dt} = -q\Omega P^{j+1}_{i} + \dfrac{\eta_T}{2} \: \left( \dfrac{Q^{j+1}_{i+1} - 2Q^{j+1}_{i} + Q^{j+1}_{i-1}}{dz^2} + \dfrac{Q^{j}_{i+1} - 2Q^{j}_{i} + Q^{j}_{i-1}}{dz^2} \right) $$

Putting $\dfrac{\alpha_0 \: dt}{2 \: dz} = \mu$ and $\dfrac{\eta_T \: dt}{2 \: dz^2} = \nu$, we separate the present time-step $(j+1)$ and the past time-step $(j)$ as $$ \left(1+2\nu \right) P^{j+1}_{i} - \nu P^{j+1}_{i+1} - \nu P^{j+1}_{i-1} + \mu Q^{j+1}_{i+1} - \mu Q^{j+1}_{i} = \left(1-2\nu \right) P^{j}_{i} + \nu P^{j}_{i+1} + \nu P^{j}_{i-1} - \mu Q^{j}_{i+1} + \mu Q^{j}_{i} $$ and $$ \left(1+2\nu \right) Q^{j+1}_{i} - \nu Q^{j+1}_{i+1} - \nu Q^{j+1}_{i-1} + dt \: q \Omega P^{j+1}_{i} = \left(1-2\nu \right) Q^{j}_{i} + \nu P^{j}_{i+1} + \nu P^{j}_{i-1} $$

Let $U = \left[ P \:\: Q\right]^T$ This whole coupled equation can be simplified to single variable $U$ as $$ \left[\begin{array}{cc}1+2 \nu & -\mu \\ dt\: q \Omega & 1+2 \nu \end{array}\right] U^{j+1}_{i} + \left[\begin{array}{cc}-\nu & \mu \\ 0 & -\nu \end{array}\right] U^{j+1}_{i+1} + \left[\begin{array}{cc}-\nu & 0 \\ 0 & -\nu \end{array}\right] U^{j+1}_{i-1} = \left[\begin{array}{cc}1-2 \nu & \mu \\ 0 & 1-2 \nu \end{array}\right] U^{j}_{i} + \left[\begin{array}{cc}\nu & -\mu \\ 0 & \nu \end{array}\right] U^{j}_{i+1} + \left[\begin{array}{cc}\nu & 0 \\ 0 & \nu \end{array}\right] U^{j}_{i-1} $$

Renaming these $2 \times 2$ matrices $E$, $F$, $G$, $H$, $I$ and $J$ from left to right respectively, we can write it as $$ E U^{j+1}_{i} + F U^{j+1}_{i+1} + G U^{j+1}_{i-1} = H U^{j}_{i} + I U^{j}_{i+1} + J U^{j}_{i-1} $$

Again we can now form the final matrix, however, this time, it is a bit different. $$ \widetilde{M}U^{j+1} = \widetilde{N}U^{j} $$ Here, the final matrices $\widetilde{M}$ and $\widetilde{N}$ are the tensor product of the above 6 matrices and the M, N matrices from the previous chapter. The final matrix would look something like this

$$ \widetilde{M} = \left[\begin{array}{cc} \left[\begin{array}{ccccc}1+2 \nu & -\nu & 0 & \cdots & 0 \\ -\nu & 1+2 \nu & -\nu & \ddots & \vdots \\ 0 & -\nu & 1+2 \nu & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & -\nu \\ 0 & \cdots & 0 & -\nu & 1+2 \nu\end{array}\right] & \left[\begin{array}{ccccc}-\mu & \mu & 0 & \cdots & 0 \\ 0 & -\mu & \mu & \ddots & \vdots \\ 0 & 0 & -\mu & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \mu \\ 0 & \cdots & 0 & 0 & -\mu \end{array}\right] \\ & \\ \left[\begin{array}{ccccc} \dfrac{dt\: q \Omega}{2} & 0 & 0 & \cdots & 0 \\ 0 & \dfrac{dt\: q \Omega}{2} & 0 & \ddots & \vdots \\ 0 & 0 & \dfrac{dt\: q \Omega}{2} & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & 0 & \dfrac{dt\: q \Omega}{2} \end{array}\right] & \left[\begin{array}{ccccc}1+2 \nu & -\nu & 0 & \cdots & 0 \\ -\nu & 1+2 \nu & -\nu & \ddots & \vdots \\ 0 & -\nu & 1+2 \nu & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & -\nu \\ 0 & \cdots & 0 & -\nu & 1+2 \nu\end{array}\right] \end{array}\right] $$

and

$$ \widetilde{N} = \left[\begin{array}{cc} \left[\begin{array}{ccccc}1-2 \nu & \nu & 0 & \cdots & 0 \\ \nu & 1-2 \nu & \nu & \ddots & \vdots \\ 0 & \nu & 1-2 \nu & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \nu \\ 0 & \cdots & 0 & \nu & 1-2 \nu\end{array}\right] & \left[\begin{array}{ccccc} \mu & -\mu & 0 & \cdots & 0 \\ 0 & \mu & -\mu & \ddots & \vdots \\ 0 & 0 & \mu & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & -\mu \\ 0 & \cdots & 0 & 0 & \mu \end{array}\right] \\ & \\ \left[\begin{array}{ccccc} -\dfrac{dt\: q \Omega}{2} & 0 & 0 & \cdots & 0 \\ 0 & -\dfrac{dt\: q \Omega}{2} & 0 & \ddots & \vdots \\ 0 & 0 & -\dfrac{dt\: q \Omega}{2} & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & 0 & -\dfrac{dt\: q \Omega}{2} \end{array}\right] & \left[\begin{array}{ccccc}1-2 \nu & \nu & 0 & \cdots & 0 \\ \nu & 1-2 \nu & \nu & \ddots & \vdots \\ 0 & \nu & 1-2 \nu & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \nu \\ 0 & \cdots & 0 & \nu & 1-2 \nu\end{array}\right] \end{array}\right] $$

The final answer i.e., $U^{j+1}$ is obtained as $$ U^{j+1} = \widetilde{M}^{-1}\widetilde{N}U^{j} $$

2.3.2 - Grid and boundary conditions

We will use the same grid structure and boumndary conditions that we used in chapter 1.

---

---

THE CODE

For the full code with al functions and plotting scripts, please visit my github repository.

Now let us get into some examples to understand how the magnetic field varies with different dynamo numbers.

from my_code import *
from plotting import *

Example 1

We have already defined all the parameters in section 2.2.2.

For this example, we define $\alpha_0 = 3 \times 10^5$. Putting this value, we get the dynamo number $D = 517.5$

def init_cond_Br(x):
    return (1-x**2)*np.cos(np.pi*x) #(1-x**2)*np.cos(2*np.pi*x)

def init_cond_Bphi(x):
    return -(1-x**2)*np.cos(np.pi*x)

def source_term(x, t):
    return 0

z = np.linspace(-1, 1, 101)
title_1 = r'Parabolic seed field '+'\n'+r'$ y = 1-x^2$'
title_2 = r'Cos square seed field '+'\n'+r'$ y = \cos^2(\pi x/2)$'
global_title = 'INITIAL CONDITIONS FOR MAGNETIC FIELD COMPONENTS'

plot_init_cond(z, init_cond_Br, init_cond_Bphi, title_1, title_2, global_title)
plt.show()

# Constants and parameters
eta_T = 3.48e-3    # magnetic diffusivity
alpha = 3e-5    # alpha effect
Omega = 40
q = -0.5
t_max = 1000     # total simulation time
z_min = -1.0     # minimum thickness of the disc
z_max = 1.0     # thickness of the disc
dt = t_max/200       # time step
dz = 0.01       # spatial step in z direction

# Spatial grid
z = np.linspace(z_min, z_max, int((z_max - z_min) / dz) + 1)
t = np.linspace(0, t_max, int(t_max / dt) + 1)

# Coefficients for the matrix A and B
rho = eta_T*dt/(2*dz**2)
sigma = alpha*dt/(2*dz)

A = matrix(len(z), 1+2*rho, -sigma, q*Omega*dt/2, 1+2*rho, -rho, sigma, 0, -rho, -rho, 0, 0, -rho)
B = matrix(len(z), 1-2*rho, sigma, -q*Omega*dt/2, 1-2*rho, rho, -sigma, 0, rho, rho, 0, 0, rho)

# Solve the diffusion equation in radial direction
solution = crank_nicolson_mod(len(z), len(t), init_cond_Br(z), init_cond_Bphi(z), A, B)

B_r = solution[:len(z), :]
B_phi = solution[len(z):, :]

# Plot the solution in imshow
plot_diff(t, z, B_r, B_phi, D=517.5)
plt.show()

Here is an animation showing the evolution of $B_r$ and $B_{\phi}$ with time.

One can see that in this example, the fields are decaying with time. We will calculate the decay rate at the midplane, i.e., $z=0$. To calculate the decay rate, we will find how the maxima of the field at the midplane varies with time.

For this, we will first calculate the local maxima of each variation and then fit an exponential decay to it. Let us first calculate the whole thing for a longer duration to get more data points.

# Constants and parameters
eta_T = 3.48e-3    # magnetic diffusivity
alpha = 3e-5    # alpha effect
Omega = 40
q = -0.5
t_max = 2000     # total simulation time
z_min = -1.0     # minimum thickness of the disc
z_max = 1.0     # thickness of the disc
dt = t_max/500       # time step
dz = 0.01       # spatial step in z direction

# Spatial grid
z = np.linspace(z_min, z_max, int((z_max - z_min) / dz) + 1)
t = np.linspace(0, t_max, int(t_max / dt) + 1)

# Coefficients for the matrix A and B
rho = eta_T*dt/(2*dz**2)
sigma = alpha*dt/(2*dz)

A = matrix(len(z), 1+2*rho, -sigma, q*Omega*dt/2, 1+2*rho, -rho, sigma, 0, -rho, -rho, 0, 0, -rho)
B = matrix(len(z), 1-2*rho, sigma, -q*Omega*dt/2, 1-2*rho, rho, -sigma, 0, rho, rho, 0, 0, rho)

# Solve the diffusion equation in radial direction
solution = crank_nicolson_mod(len(z), len(t), init_cond_Br(z), init_cond_Bphi(z), A, B)

B_r = solution[:len(z), :]
B_phi = solution[len(z):, :]

B_mid_r = B_r[int(len(z)/2), :]
B_mid_phi = B_phi[int(len(z)/2), :]

B_mid = np.sqrt(B_mid_r**2 + B_mid_phi**2)

plt.figure(figsize=(6, 4))
plt.plot(t, B_mid)
x_max, y_max = find_local_maxima(t, B_mid)
plt.plot(x_max, y_max, 'ro')
plt.xlabel('Time')
plt.ylabel('Magnetic field')
plt.title('Magnetic field at the midplane')
plt.show()

decay_rate = get_decay_rate(t, B_mid)
if decay_rate < 0:
    print(r"decay rate at the mid plane:", format(-decay_rate, ".3e"))
else:
    print(r"growth rate at the mid plane:", format(decay_rate, ".3e"))

print('Dynamo number ', alpha*1.725e7)

decay rate at the mid plane: 1.607e-04
Dynamo number  517.5

Example 2

For this example, we choose $\alpha_0 = 4 \times 10^5$. Putting this value, we get the dynamo number $D = 690$

def init_cond_Br(x):
    return (1-x**2)*np.cos(np.pi*x) #(1-x**2)*np.cos(2*np.pi*x)

def init_cond_Bphi(x):
    return -(1-x**2)*np.cos(np.pi*x)

def source_term(x, t):
    return 0

z = np.linspace(-1, 1, 101)
title_1 = r'Parabolic seed field '+'\n'+r'$ y = 1-x^2$'
title_2 = r'Cos square seed field '+'\n'+r'$ y = \cos^2(\pi x/2)$'
global_title = 'INITIAL CONDITIONS FOR MAGNETIC FIELD COMPONENTS'

plot_init_cond(z, init_cond_Br, init_cond_Bphi, title_1, title_2, global_title)
plt.show()

# Constants and parameters
eta_T = 3.48e-3    # magnetic diffusivity
alpha = 4e-5    # alpha effect
Omega = 40
q = -0.5
t_max = 1000     # total simulation time
z_min = -1.0     # minimum thickness of the disc
z_max = 1.0     # thickness of the disc
dt = t_max/200       # time step
dz = 0.01       # spatial step in z direction

# Spatial grid
z = np.linspace(z_min, z_max, int((z_max - z_min) / dz) + 1)
t = np.linspace(0, t_max, int(t_max / dt) + 1)

# Coefficients for the matrix A and B
rho = eta_T*dt/(2*dz**2)
sigma = alpha*dt/(2*dz)

A = matrix(len(z), 1+2*rho, -sigma, q*Omega*dt/2, 1+2*rho, -rho, sigma, 0, -rho, -rho, 0, 0, -rho)
B = matrix(len(z), 1-2*rho, sigma, -q*Omega*dt/2, 1-2*rho, rho, -sigma, 0, rho, rho, 0, 0, rho)

# Solve the diffusion equation in radial direction
solution = crank_nicolson_mod(len(z), len(t), init_cond_Br(z), init_cond_Bphi(z), A, B)

B_r = solution[:len(z), :]
B_phi = solution[len(z):, :]

# Plot the solution in imshow
plot_diff(t, z, B_r, B_phi, D=690)
plt.show()

Here is an animation showing the evolution of $B_r$ and $B_{\phi}$ with time.

One can see that in this example, the fields are growing with time. We will calculate the growth rate at the midplane, i.e., $z=0$. To calculate the growth rate, we will find how the maxima of the field at the midplane varies with time.

For this, we will first calculate the local maxima of each variation and then fit an exponential growth to it. Let us first calculate the whole thing for a longer duration to get more data points.

# Constants and parameters
eta_T = 3.48e-3    # magnetic diffusivity
alpha = 4e-5
Omega = 40
q = -0.5
t_max = 2000     # total simulation time
z_min = -1.0     # minimum thickness of the disc
z_max = 1.0     # thickness of the disc
dt = t_max/500       # time step
dz = 0.01       # spatial step in z direction

# Spatial grid
z = np.linspace(z_min, z_max, int((z_max - z_min) / dz) + 1)
t = np.linspace(0, t_max, int(t_max / dt) + 1)

# Coefficients for the matrix A and B
rho = eta_T*dt/(2*dz**2)
sigma = alpha*dt/(2*dz)

A = matrix(len(z), 1+2*rho, -sigma, q*Omega*dt/2, 1+2*rho, -rho, sigma, 0, -rho, -rho, 0, 0, -rho)
B = matrix(len(z), 1-2*rho, sigma, -q*Omega*dt/2, 1-2*rho, rho, -sigma, 0, rho, rho, 0, 0, rho)

# Solve the diffusion equation in radial direction
solution = crank_nicolson_mod(len(z), len(t), init_cond_Br(z), init_cond_Bphi(z), A, B)

B_r = solution[:len(z), :]
B_phi = solution[len(z):, :]

B_mid_r = B_r[int(len(z)/2), :]
B_mid_phi = B_phi[int(len(z)/2), :]

B_mid = np.sqrt(B_mid_r**2 + B_mid_phi**2)

x_max, y_max = find_local_maxima(t, B_mid)
plot_aO_decay(t, B_mid, x_max, y_max)
plt.show()

decay_rate = get_decay_rate(t, B_mid)
if decay_rate < 0:
    print(r"decay rate at the mid plane:", format(-decay_rate, ".3e"))
else:
    print(r"growth rate at the mid plane:", format(decay_rate, ".3e"))

print('Dynamo number ', alpha*1.725e7)

growth rate at the mid plane: 2.167e-03
Dynamo number  690.0

Calculation of Critical Dynamo number $D_c$

# Constants and parameters
eta_T = 3.48e-3    # magnetic diffusivity
Omega = 40
q = -0.5
t_max = 2000     # total simulation time
z_min = -1.0     # minimum thickness of the disc
z_max = 1.0     # thickness of the disc
dt = t_max/500       # time step
dz = 0.01       # spatial step in z direction

tol = 1e-6

def f(D):
    alpha = D/1.725e7
    # Spatial grid
    z = np.linspace(z_min, z_max, int((z_max - z_min) / dz) + 1)
    t = np.linspace(0, t_max, int(t_max / dt) + 1)

    # Coefficients for the matrix A and B
    rho = eta_T*dt/(2*dz**2)
    sigma = alpha*dt/(2*dz)

    A = matrix(len(z), 1+2*rho, -sigma, q*Omega*dt/2, 1+2*rho, -rho, sigma, 0, -rho, -rho, 0, 0, -rho)
    B = matrix(len(z), 1-2*rho, sigma, -q*Omega*dt/2, 1-2*rho, rho, -sigma, 0, rho, rho, 0, 0, rho)

    # Solve the diffusion equation in radial direction
    solution = crank_nicolson_mod(len(z), len(t), init_cond_Br(z), init_cond_Bphi(z), A, B)

    B_r = solution[:len(z), :]
    B_phi = solution[len(z):, :]
    B_mid_r = B_r[int(len(z)/2), :]
    B_mid_phi = B_phi[int(len(z)/2), :]

    B_mid = np.sqrt(B_mid_r**2 + B_mid_phi**2)

    decay_rate = get_decay_rate(t, B_mid)

    return decay_rate

print('Calculating critical Dynamo number...')
D_c = bisection(f, 200, 800, tol)
print('Critical Dynamo number Dc = ', D_c)

Calculating critical Dynamo number...
Critical Dynamo number Dc =  528.857421875